include 'mkl_vsl.f90'

module modising2d
    use modinip
    use ifport
    USE MKL_VSL_TYPE
    USE MKL_VSL
    implicit none
    private

    ! ------------------------------------------------------
    ! Struktura z parametrami ukladu
    ! ------------------------------------------------------
    type global_parameters
        double precision :: T   ! temperatura
        double precision :: B   ! pole magnetyczne
        double precision :: Jxy ! oddzialywanie w plaszczyznie
        double precision :: Jz  ! oddzialywanie w pionie
        integer          :: nx  ! rozmiar siatki
    end type global_parameters
    ! ------------------------------------------------------
    ! Struktura z danymi siatki
    ! ------------------------------------------------------
    type lattice_parameters
        double precision :: energy
        double precision :: mean_energy     ! <E>
        double precision :: mean_sqr_energy ! <E^2>
        double precision :: Cw              ! Cieplo wlasciwe: <E^2> - <E>^2
        double precision :: magnetization          ! M   = <S>/N   - sredni spin na wezel
        double precision :: mean_magnetization     !  <M>
        double precision :: mean_sqr_magnetization !  <M^2>
        double precision :: Xi              ! Podatnosc: <M^2> - <M>^2

        integer          :: no_rejections ! liczba odrzuconych krokow Metropolisa
        doubleprecision  :: acc_ratio     ! stosunek zaakceptowanych krokow do calkowitej liczby
    end type lattice_parameters

    ! ------------------------------------------------------
    !
    ! ------------------------------------------------------
    type(global_parameters)                    :: g_params
    type(lattice_parameters)                   :: l_params
    integer,dimension(:),allocatable           :: lattice
    doubleprecision,dimension(:),allocatable   :: sublattice ! vektor ze spinami substratu
    doubleprecision,dimension(:),allocatable   :: avelattice ! vektor ze srednimi spinami na sieci
    doubleprecision,dimension(:),allocatable   :: randvec    ! vektor z losowymi liczbami o rozmiarze lattice
    integer,dimension(:),allocatable           :: randpos    ! vektor z losowymi indeksami przy iteracji MC
    integer,dimension(:),allocatable   	       :: border     ! vektor - "stos" z indeksami spinow odwroconych do algorytmu Wolffa
    integer,dimension(:),allocatable   	       :: cluster  	!tablica w ktorej sa indeksy atomow klastra
								!sprawdzamy zmiane energii deltaE oddz. z dolna siecia
								!i odwracamy podobnie jak w metropolisie - jesli rand()<exp(-deltaE/T)
    ENUM,BIND(C)
          ENUMERATOR :: ZAPISZ_SIEC         = 1
          ENUMERATOR :: ZAPISZ_PODSIEC      = 2
          ENUMERATOR :: ZAPISZ_SREDNIA_SIEC = 3
    END ENUM
    ! zmienne zwiazane z generowaniem liczb losowych
    TYPE (VSL_STREAM_STATE)                    :: stream
    integer,parameter :: brng   = VSL_BRNG_MT19937
    integer,parameter :: method = VSL_RNG_METHOD_UNIFORM_STD
    integer(kind=4)   :: errcode ! wyscie bledu dla mkl vsl
    ! ----------------------------------------------------
    public :: g_params , l_params
    public :: ising_wczytaj_parametry , ising_alokuj_pamiec , ising_dealokuj_pamiec , ising_inicjalizuj_siec
    public :: ising_zapisz_do_pliku , ising_zapisz_przekroj_X, ising_metropolis_cycle , ising_cycle
    public :: lattice
    public :: ZAPISZ_SIEC , ZAPISZ_PODSIEC , ZAPISZ_SREDNIA_SIEC
    contains

    ! ----------------------------------------------
    ! Wczytujemy ustawienia z pliku config.ini
    ! tj globalne ustawienia symulacji takie jak
    ! temperatura, pole B, czy np. tryb symulacji
    ! ----------------------------------------------
    subroutine ising_wczytaj_parametry(filename)
    character(*), intent(in) :: filename
        print*,"-----------------------------------------"
        print*,"   Wczytywanie danych z config.ini:"
        print*,"-----------------------------------------"
        call setIniFilename(filename)
        call getDoubleValue("Dane","T",g_params%T)
        call getDoubleValue("Dane","B",g_params%B)
        call getDoubleValue("Dane","Jxy",g_params%Jxy)
        call getDoubleValue("Dane","Jz",g_params%Jz)
        call getIntValue   ("Dane","nx",g_params%nx)
    end subroutine ising_wczytaj_parametry

    ! -----------------------------------------------
    ! Alokacja pamieci tablic
    ! -----------------------------------------------
    subroutine ising_alokuj_pamiec()
        integer :: m
        integer :: seed(100)
        if(allocated(lattice)) then
            deallocate(lattice)
        endif
        if(allocated(randvec)) then
            deallocate(randvec)
        endif
        if(allocated(randpos)) then
            deallocate(randpos)
        endif
        if(allocated(sublattice)) then
            deallocate(sublattice)
        endif
        if(allocated(avelattice)) then
            deallocate(avelattice)
        endif
        if(allocated(border)) then
            deallocate(border)
        endif
        if(allocated(cluster)) then
            deallocate(cluster)
        endif
        allocate(lattice   (g_params%nx*g_params%nx))
        allocate(randvec   (g_params%nx*g_params%nx))
        allocate(randpos   (g_params%nx*g_params%nx))
        allocate(sublattice(g_params%nx*g_params%nx))
        allocate(avelattice(g_params%nx*g_params%nx))
        allocate(border    (g_params%nx*g_params%nx))
        allocate(cluster   (g_params%nx*g_params%nx))

        CALL RANDOM_SEED
        CALL RANDOM_SEED(SIZE=M)
        CALL RANDOM_SEED(GET=seed(1:M))
        print*,"Poczatkowy SEED:",sum(seed(:))
        errcode = vslnewstream( stream, brng,  sum(seed(:)) )

        if(errcode  /= 0 ) then
            print*,"MKL: Blad inicjalizacji VSL - ERRCODE:",errcode
        endif
        avelattice = 0
    end subroutine ising_alokuj_pamiec

    subroutine ising_dealokuj_pamiec()
        if(allocated(lattice))    deallocate(lattice)
        if(allocated(randvec))    deallocate(randvec)
        if(allocated(randpos))    deallocate(randpos)
        if(allocated(sublattice)) deallocate(sublattice)
        if(allocated(avelattice)) deallocate(avelattice)
        if(allocated(avelattice)) deallocate(border)
        if(allocated(avelattice)) deallocate(cluster)
        errcode=vsldeletestream( stream )
    end subroutine ising_dealokuj_pamiec

    ! -------------------------------------------------------------
    ! Inicjalizacja sieci losowymi liczbami -1,+1 liczbami.
    ! oder_param = 0.0 - cala siec -1
    ! oder_param = 0.5 - calkowity nieporzadek
    ! oder_param = 1.0 - cala siec +1
    ! -------------------------------------------------------------
    subroutine ising_inicjalizuj_siec(order_param)
        double precision :: order_param
        integer          :: i , j, nx , k , kpow
        call ising_wypelnij_randvec(0.0D0,1.0D0)
        nx = g_params%nx
        lattice = -1
        do i = 1 , nx*nx
            if( randvec( i ) < order_param ) lattice(i) = 1
        enddo
        l_params%energy = ising_calc_energy()
        l_params%mean_energy      = 0
        l_params%mean_sqr_energy  = 0
        sublattice = 1
        kpow = 1
        do i = 1 , nx
        do j = 1 , nx
            if(  abs(j-nx/2.0)**kpow < (nx/4.0)**kpow ) sublattice(to2D(i,j)) = -1
        enddo
        enddo

    end subroutine ising_inicjalizuj_siec



    ! --------------------------------------------------------
    ! Funkcja pomocnicza sluzaca do generowania ciagu liczb
    ! losowych z przedzialu od xmin do xmax. Wypelniana jest
    ! cala tablica randvec. Ma ona tyle elementow co siec isinga.
    ! --------------------------------------------------------
    subroutine ising_wypelnij_randvec(xmin,xmax)
        doubleprecision :: xmin, xmax
        integer :: errcode
        errcode = vdrnguniform( method, stream, size(randvec), randvec, xmin, xmax);
    end subroutine ising_wypelnij_randvec
    subroutine ising_wypelnij_randpos()
        integer :: errcode
        errcode = virnguniform( method, stream, size(randpos), randpos, 1, size(randpos)+1);
    end subroutine ising_wypelnij_randpos

    ! --------------------------------------------------------
    ! Zapisuje aktualny stan sieci to pliku:
    ! do wyboru mamy rozne sieci do zapisania:
    !    siec         - ZAPISZ_SIEC,
    !    podsiec      - ZAPISZ_PODSIEC
    !    srednia siec - ZAPISZ_SREDNIA_SIEC
    ! --------------------------------------------------------
    subroutine ising_zapisz_do_pliku(filename,enum_zapisz)
        character(*),intent(in) :: filename
        integer,intent(in)      :: enum_zapisz
        integer :: i,nx
        nx = g_params%nx
        open(unit=111,file=filename)
        write(111,*),nx
        do i = 1 , nx*nx
            select case(enum_zapisz)
                case (ZAPISZ_SIEC)
                    write(111,*),lattice(i)
                case (ZAPISZ_PODSIEC)
                    write(111,*),sublattice(i)
                case (ZAPISZ_SREDNIA_SIEC)
                    write(111,*),avelattice(i)
            endselect
        enddo
        close(111)
    end subroutine ising_zapisz_do_pliku

    ! -------------------------------------------------------
    ! Zapisuje do pliku jedna linie dla zadanej pozycji Y.
    ! -------------------------------------------------------
    subroutine ising_zapisz_przekroj_X(filename,ypos)
        character(*),intent(in) :: filename
        integer,intent(in)      :: ypos
        doubleprecision         :: ysum1 , ysum2 , ysum3
        integer :: i,j
        open(unit=111,file=filename)
        do i = 2 , g_params%nx-1
            ysum1 = 0
            ysum2 = 0
            ysum3 = 0
            do j = 1 , g_params%nx
                ysum1 = ysum1 + (avelattice(to2D(j,i)))
                ysum2 = ysum2 + (avelattice(to2D(j,i+1)))
                ysum3 = ysum3 + (avelattice(to2D(j,i-1)))
            enddo
            write(111,"(i10,20e20.8)")i,dble(lattice(to2D(ypos,i))),&
            & dble(sublattice(to2D(ypos,i))),avelattice(to2D(ypos,i)),ysum1/g_params%nx,&
            & abs(ysum2 - ysum3)/(2.0*g_params%nx)
        enddo
        close(111)
    endsubroutine ising_zapisz_przekroj_X

    ! ----------------------------------------------------------------
    ! Wykonuje jeden cykl metropolisa - to jest tyle zamian spinow
    ! co liczba spinow w siatce.
    ! ----------------------------------------------------------------
    subroutine ising_metropolis_cycle()
        integer :: jpos , iter , nx , L , i
        integer :: siteL , siteR , siteU , siteD
        double precision :: deltaE , invT

        call ising_wypelnij_randvec(0.0D0,1.0D0) ! generacja wektora losowych liczb
        call ising_wypelnij_randpos()            ! generacja wektora polozen losowych

        nx   = g_params%nx
        L    = nx*nx
        invT = 1.0/g_params%T

        do iter = 1 , L
            i     = randpos(iter)
            jpos  = ( ( i - 1 ) / NX )*NX + 1

            siteL = mod( i + NX     , NX) + jpos
            siteR = mod( i - 2 + NX , NX) + jpos
            siteU = mod( (i - 1 + nx + L) , L) + 1
            siteD = mod( (i - 1 - nx + L) , L) + 1

            deltaE = 2 * g_params%Jxy * lattice(i)*(lattice(siteU)+lattice(siteD)+lattice(siteL)+lattice(siteR))
            deltaE = deltaE + 2*g_params%Jz*lattice(i)*sublattice(i) !czy tu nie ma byc *2?

            if( randvec(iter) < exp(-deltaE * invT) ) then
                lattice(i) = -lattice(i)
            else
                l_params%no_rejections = l_params%no_rejections + 1
            endif
        enddo ! end of iter = 1 , L

        avelattice                = avelattice + lattice
        ! energia
        l_params%energy           = ising_calc_energy()
        l_params%mean_energy      = l_params%mean_energy     + l_params%energy
        l_params%mean_sqr_energy  = l_params%mean_sqr_energy + l_params%energy**2
        ! magnetyzacja
        l_params%magnetization          = dble(sum(lattice))/dble(L)
        l_params%mean_magnetization     = l_params%mean_magnetization     + abs(l_params%magnetization)
        l_params%mean_sqr_magnetization = l_params%mean_sqr_magnetization + l_params%magnetization**2


    end subroutine ising_metropolis_cycle


    ! ----------------------------------------------------------------
    ! Wykonuje jeden cykl wolffa - zamiana spinow w klastrach
    ! ----------------------------------------------------------------
    subroutine ising_wolff_cycle(sum_cluster)
    	integer              ::  sum_cluster  !liczba atomow w klastrze
	integer 	     ::  bindex,i,jpos,nx !iter,
	double precision     ::  invT
	integer 	     ::  oldspin,newspin
	integer		     ::  site, siteL, siteR, siteU, siteD
	double precision     ::  Padd, deltaE
	double precision     ::  L
	!integer 	     :: sum_checked  !liczba spinow sprawdzonych przez algorytm - do
				!zliczania operacji aby zlozonosc Wolffa i Metropolisa byla podobna

	nx = g_params%nx
    invT = 1.0/g_params%T
	Padd = ( 1 - exp(-2.*invT)  )
	bindex = 0
	L=nx*nx
	!sum_checked=0

	    i = L*rand()+1 !indeks wylosowanego spinu
	    site = i

	    bindex = 1
	    border(bindex) = i
	    oldspin =  lattice(i)
	    newspin = -lattice(i)
	    lattice(i) = newspin
	    sum_cluster=0

		    do while( bindex .ge. 1 )
			site   = border(bindex)
			bindex=bindex-1
			sum_cluster=sum_cluster+1
			cluster(sum_cluster)=site

			jpos = ( ( site - 1 ) / NX )*NX + 1
		    	siteL = mod( site - 2 + NX , NX) + jpos
			if ( lattice(siteL) == oldspin ) then
			    if( rand()<Padd ) then
				bindex=bindex+1
				border(bindex) = siteL
				lattice(siteL) = newspin
			    endif
			endif

		    	siteR = mod( site + NX     , NX) + jpos
			if ( lattice(siteR) == oldspin ) then
			    if( rand()<Padd ) then
				bindex=bindex+1
				border(bindex) = siteR
				lattice(siteR) = newspin
			    endif
			endif

		    	siteU = mod( (site - 1 + nx + L) , L) + 1
			if ( lattice(siteU) == oldspin ) then
			    if( rand()<Padd ) then
				bindex=bindex+1
				border(bindex) = siteU
				lattice(siteU) = newspin
			    endif
			endif

		    	siteD = mod( (site - 1 - nx + L) , L) + 1
			if ( lattice(siteD) == oldspin ) then
			    if( rand()<Padd ) then
				bindex=bindex+1
				border(bindex) = siteD
				lattice(siteD) = newspin
			    endif
			endif

		    enddo! end of while bindex .ge. 1


        avelattice = avelattice   + lattice
        l_params%energy           = ising_calc_energy()
        l_params%mean_energy      = l_params%mean_energy     + l_params%energy
        l_params%mean_sqr_energy  = l_params%mean_sqr_energy + l_params%energy**2
        ! magnetyzacja
        l_params%magnetization          = dble(sum(lattice))/dble(L)
        l_params%mean_magnetization     = l_params%mean_magnetization     + abs(l_params%magnetization)
        l_params%mean_sqr_magnetization = l_params%mean_sqr_magnetization + l_params%magnetization**2

    end subroutine ising_wolff_cycle

    ! ---------------------------------------------------------
    ! Wykonuje zadana liczbe cykli -no_iter metropolisja.
    ! Cykl poprzedzony jest nie zliczana liczba cykli termalizacji
    ! ukladu.
    ! ---------------------------------------------------------
    subroutine ising_cycle(no_iter,therm_time)
        integer :: no_iter,therm_time,sum_cluster,c_iter,w_iter
        integer :: i
	integer :: method !0-metropolis, 1- wolff
        ! termalizacja
        do i = 1, therm_time
            call ising_metropolis_cycle()
        enddo

        ! zerowanie srednich wartosci
        l_params%energy                 = 0
        l_params%mean_energy            = 0
        l_params%mean_sqr_energy        = 0
        l_params%no_rejections          = 0
        l_params%magnetization          = 0
        l_params%mean_magnetization     = 0
        l_params%mean_sqr_magnetization = 0
        l_params%Xi                     = 0
        avelattice                      = 0
        ! wlasciwa symulacja
        w_iter = 0
        do i = 1, no_iter

	method=1
	if(method==0) then
           w_iter = w_iter + 1
           call ising_metropolis_cycle()
	else
            c_iter = 0
            do while(c_iter < g_params%nx**2/2)
                call ising_wolff_cycle(sum_cluster) !todo: wybor wolffa lub metr.
                c_iter = c_iter + sum_cluster
                w_iter = w_iter + 1
            enddo
	endif

		!print*,l_params%mean_magnetization, l_params%magnetization, w_iter

        enddo
        ! obliczamy srednie wartosci

        avelattice               =  avelattice / w_iter
        l_params%mean_energy     =  l_params%mean_energy    / w_iter
        l_params%mean_sqr_energy =  l_params%mean_sqr_energy/ w_iter
        l_params%Cw              = (l_params%mean_sqr_energy - l_params%mean_energy**2)/(g_params%T)**2/(g_params%nx)**2
		!print*,""!l_params%mean_magnetization, l_params%magnetization, w_iter
        ! magnetyzacja i podatnosc
        l_params%mean_magnetization     =  l_params%mean_magnetization    /w_iter
        l_params%mean_sqr_magnetization =  l_params%mean_sqr_magnetization/w_iter
        l_params%Xi                     = (l_params%mean_sqr_magnetization - l_params%mean_magnetization**2)/(g_params%T)

        l_params%acc_ratio       = ( 1 - dble(l_params%no_rejections)/dble(w_iter*g_params%nx**2) )
    end subroutine ising_cycle


    ! ----------------------------------------------
    ! Dokonuje mapowania z ukladu numerowanego dwoma
    ! indeksami (i,j) do ukladu numerowanego jednym
    ! indeksem i.
    ! ----------------------------------------------
    integer function to2D(i,j)
        integer :: i,j
        to2D = i + (j-1)*g_params%nx
    end function to2D

    ! --------------------------------------------------------
    ! Oblicza energie ukladu
    ! --------------------------------------------------------
    doubleprecision function ising_calc_energy()
        doubleprecision :: E
        integer :: i,j,nx
        E = 0
        nx = g_params%nx
        do i = 1 , nx
        do j = 1 , nx
            E = E - g_params%Jxy*lattice(to2D(i,j))*( lattice(To2D(mod(i,nx)+1,j)) + lattice(To2D(i,mod(j,nx)+1)) )
            E = E + g_params%B*lattice(To2D(i,j))
            E = E - g_params%Jz*lattice(to2D(i,j))*sublattice(to2D(i,j))
        enddo
        enddo
        ising_calc_energy = E
    end function ising_calc_energy





end module modising2d

